%% 清空环境
clear; close all; clc;

%% 参数设置
a0 = -0.82;      % 不同涡旋对应不同 a0
a1 = 0.83;
a2 = 0.14;
b0 = 1.3;
b1 = 0.1;
b2 = 1;
d1 = 0;
d2 = 0;

%% 初始条件
x0 = 0.1;
y0 = 0.1;
z0 = 0.1;

%% 迭代总数
N = 500000;
burn_in = 1000;

%% 预分配数组
x = zeros(1, N);
y = zeros(1, N);
z = zeros(1, N);

% 初值
x(1) = x0;
y(1) = y0;
z(1) = z0;

%% 开始迭代
for n = 1:N-1
    x(n+1) = y(n) + x(n);
    y(n+1) = sin( y(n)*a0*sin(a1 + a2*x(n))*y(n) ) + b2*sin(z(n)) + d1;
    z(n+1) = b0*y(n) + b1*z(n) + d2;
end

%% 绘图 (3D 相图)
figure;
plot3( x(burn_in:end), y(burn_in:end), z(burn_in:end), '.', 'MarkerSize', 1 );
xlabel('x_n');
ylabel('y_n');
zlabel('z_n');
title('Fig.6 三维相图');
grid on;
view(-45,30);  % 调整视角